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A system of active colloidal particles driven by harmonic potentials to oscillate about the vertices 
of a regular polygon, with hydrodynamic coupling between all particles, is described by a piece-wise 
linear model which exhibits various patterns of synchronization. Analytical solutions are obtained 
for this class of dynamical systems. Depending only on the number of particles, the synchronization 
occurs into states in which nearest neighbors oscillate either in-phase, or anti-phase, or in phase- 
locked (time-shifted) trajectories. 



I. INTRODUCTION 

In passive colloidal particle systems, hydrodynamic in- 
teractions determine the dissipative and dynamical be- 
havior, rather than the real-space conformation; as such 
they influence the response properties rather than the 
equilibrium behavior. However in actively driven sys- 
tems that settle into a steady state, it can be the hydro- 
dynamic interactions that determine the characteristics 
of the steady state. For example, in biological flows actu- 
ated by cilia [T], the hydrodynamic coupling is possibly 
the most important feature involved in the synchroniza- 
tion of cilia beats [2] . 

In this paper we analyze theoretically a very simple dy- 
namical system which was recently proven to adequately 
describe the motion of a system of colloidal spheres, 
maintained in oscillations by optical traps (31 0]. Mo- 
tion occurs at low Reynolds number [5], and hydrody- 
namic forces between particles depend on the relative 
velocity Oil]. 

The remarkable feature of the model is that the hy- 
drodynamic interaction of the spheres leads to the syn- 
chronization of their phases, for any number of beads, 
for all initial conditions and for generic values of physical 
parameters (sphere size, viscosity of the fluid, stiffness 
of the forces). The synchronized state depends on the 
number of spheres, and on wether the optical traps pro- 
vide attractive or repulsive forces. This produces a rich 
pattern of synchronized configurations. 

The system of equations of motion is piece- wise linear. 
It has periodic solutions, which are linear combinations 
of the normal modes. Whilst we do not have a mathe- 
matical proof which would explain why this dynamical 
system always converges to a synchronized solution, we 
have performed a normal mode analysis and we show 
that based on this it is possible to predict which synchro- 
nized solution is dominant. This analytical argument 
is in agreement with extensive numerical integration of 
the deterministic (no thermal noise) system, and with 
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the experiments and simulations reported in the com- 
panion paper which include the effect of thermal noise [1] . 

Let us outline the physical system and its model, 
leading to the dynamic deterministic system of equation 
of motion which will be our main concern in the next 
section. We refer the reader to the literature which 
details a series of experiments and their analysis. By 
focused laser beams, it was possible to confine effectively 
an array of micron-sized beads, each one into a harmonic 
well. These colloidal particles interact exclusively 
through the hydrodynamic coupling. The collective fluc- 
tuation modes resulting from hydrodynamic interactions 
have been studied in linear arrays [8| and ring arrays 
[HI Hn] • The "geometric switch" active driving model we 
study here was proposed in the context of cilia by ilj; it 
was studied analytically for two beads and for infinite 
linear chains by Cosentino Lagomarsino and Bassetti 
in [TT], where they showed that an infinite linear chain 
of oscillators can sustain a traveling wave solution. This 
model was then realized experimentally in [3], by using 
optical tweezers to maintain a pair of colloidal spheres 
in oscillation by switching the positions of optical traps 
when a sphere reaches its limit position. This rule leads 
to oscillations that are bounded in amplitude but free in 
phase and period. The interaction between the oscilla- 
tors is only through the hydrodynamic flow induced by 
their motion; the colloidal particles in the experiment 
are subject to Brownian thermal fluctuations as would 
also be found in a biological context. This system of 
two spheres in harmonic traps leads to synchronized 
motion in anti-phase |3, . A companion paper to this one 
reports experimental results and the effect of noise, in 
the system of actively moving traps for many particles 
4J. Very recently, WoUin and Stark [T2] re-considered 
this model for the case of a chain of driven oscillators, 
performing numerical simulations showing that for 
various potentials and forms of coupling the dynamic 
system evolves to a synchronized state which may be 
interpreted as a (discretized) wave propagating through 
the chain. The form of the wave depends on the driving 
forces and their being attractive or repulsive. This 
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chain is closely related to the array of driven oscillators 
studied in this paper and experimentally observed [4]. 
The main differences are the geometry (our oscillators 
are tangential to a ring) and the symmetry (in our 
model, the driving force acts in the same way in the 
two phases of each cycle, i.e. as each oscillator moves 
forward or backward). 

The periodic solutions analytically obtained in this 
paper may be interpreted as a wave propagating along 
the border of the ring, clockwise or anti-clockwise. 
Indeed if Xj{t) describes the motion of the j*''-oscillator, 
our periodic solutions have the form Xj+i(t + A) = Xj(t), 
where A is a fixed, positive or negative, time interval, 
independent on j. In a recent paper |13j Uchida and 
Golestanian derived generic conditions on the driven 
oscillators to synchronize due to the hydrodynamic 
interactions. However, as it was already noted in |12j . 
it is not clear if such conditions would apply to the 
"geometric switch" model. 

The paper is organized as follows: Section [Til outlines 
the model and the associated dynamical system. Sec- 
tion |III| shows the normal mode analysis that allows 
(in principle) a completely analytical determination of 
the periodic solutions of the system for any number of 
particles. 



Section IV describes the periodic solutions, 
if the number of particles is 3, 4, 5 and the dominant 
solution for any number of particles. In Section [V] we 
replace the attractive forces with repulsive ones and 
obtain periodic solutions trivially related to the previous 
ones. 



II. THE MODEL 

Let us first describe the system in absence of hydro- 
dynamic coupling between particles. Fig. [l|a) shows a 
ring of radius R, with n spheres depicted at the ver- 
tices of a regular n— polygon. Each vertex is the cen- 
ter of small oscillations of a particle bound to the ring. 
The j'^'-particle performs oscillations about the position 

Rj = ROf = R^U - 1), with j = 1,2,... ,n. On 
each side of this position, the amplitude of the oscilla- 
tion is A/2 — ^. This is assumed to be so small that 
we may replace the path on the arc with a segment of 
straight line, tangent to the ring, and then consider the 
motion of the j*''-particle on the tangent with a variable 
Xj: —(A/2 — ^) < Xj{t) < A/2 — ^. The oscillations are 
driven by an attractive harmonic potential (or repulsive 
harmonic potential, as discussed later in section |v| with 
center at s = ±A/2. In the hydrodynamic regime where 
Reynolds number is much smaller then unity, the inertial 
term may be neglected. 

An oscillation is set up that has fixed amplitude, but 
free phase and period. This is realized through a "ge- 
ometric switch" condition. First the trap is set at A/2, 
and particle moves along the positive x axis according to 





FIG. 1. (Color online) (a) Diagram of the ring arrangement 
with n = 6 driven oscillators, illustrating the notation used in 
the main text. Colloidal particles move tangential to the ring, 
with the center of oscillations being the vertices of a regular 
n— polygon. The local tangential co-ordinate system is shown 
for one of the top particles. The particle configuration (solid 
disks) corresponds to an in-phase condition, when all 6 par- 
ticles have just reached the geometric switch boundary. The 
particles are driven by harmonic traps towards the minimum 
of the potential, but they will reach first the other geomet- 
ric switch boundary, when at the position indicated with open 
disks. This work also considers harmonic potentials with neg- 
ative stiffness (i.e. repulsive) in section |v] (b) Trajectory of 
a driven oscillator with attractive driving force, uncoupled to 
other oscillators; x{t) is measured along the tangential direc- 
tion. 



the differential equation 



7oi(t) + K ( xit) - 2 ) - = 0' 



where k is the stiffness of the harmonic force, tq — ^q/k 
is the relaxation time, and / is a stochastic force due to 
thermal agitation of the fluid, which will be neglected in 
this work. As the particle reaches the position f — the 
attractive harmonic potential is moved at the position 
— ^. The particle inverts its overdamped motion until it 
reaches the position — ^ + ^. At this time, the poten- 
tial jumps again. The resulting oscillatory motion (this 
is either a single bead, or a bead in a system with no 
coupling), neglecting any stochastic force, is depicted in 
Fig. . The period Tq of these uncoupled oscillations 
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In a system of n beads, at any given time, the set of 
positions of the particles is a n— th dimensional vector 
x{t) = {xi{t),X2{t), ■ ■ ■ ,Xn{t)}. In the absence of the 
hydrodynamic coupling, each particle would perform the 
same periodic motion described above, a set of periodic 
overdamped relaxations. For any pair of particles, 
the "phase differences" Xj{t) — Xk{t) would be set by 
initial conditions (typically some random values) and be 
constant in time. 

Let us now include the hydrodynamic interaction. 
The motion of the j'*''-particle originates a force on the 



■th 



particle /^^ = (iJ-^) 



, where H is the Oseen 



tensor [S] . In agreement with previous analysis [SHTU] we 
are led to the system of equations 



-1 dfj(t) 
j dt 



+m = 

0, t{e) - 



0, i = l,2- 
- sin 
cos 9 



(1) 



The deterministic force Fi acting on the «-th particle 
is a harmonic force, tangent to the ring, with the "geo- 
metric switch" rule. 



x,{t) T 



A 



27r(i - 1) 



where ± f is the coordinate of the bottom of the harmonic 
well. 

The stochastic force fi{t) in Eq. [T] represents the ther- 
mal noise on the z-th particle, and it can be assumed that 

< Mt) >= 0, < Mh)fjit2) >= 2(/?)-i (if-i).j <5(ti - 

^2), f3^^ — ksT 8 . t{fi) is a versor tangent to the ring, 
at the position fi , with anti-clockwise direction. For 
small oscillations, the Oseen tensor can be approximated 
by inserting the fixed distances among centers of oscilla- 
tions, then: 



10 H"^ = SrjSafi + (1 - %) 



where ( ■'^ 
'''jy 



3a 

4r„ 



R 



r-r. ■ 



cos(j — l)27r/?i 
sin(j — l)27r/n 



with ~Rf~Rf^ ~r^,, and a = 1, 2. 

We project the equations of the linearized Langevin- 
Oseen system along the tangents, replace the global 
coordinates of the j*''-particle f^-'^ with the local one- 
dimensional coordinate x^^^ and neglect the stochastic 
force, to obtain the deterministic linear system which 
holds in the time interval where no particle hits its left 
or right boundary ±(| ~ 0'- 



Initial conditions 
x(to) 



Normal mode solution Updated initial condition 

propagated until first switch and switch state 
x(t);to<t<ti x(ti) 

s(0) s(1) 



^ Study the iterative discrete map {tj, x(ti), s^')} i=i,2... 

FIG. 2. (Color online) A solution of the dynamical system can 
be obtained by iterating the linear evolution of the system, 
which is simple in each of the eigenmodes, up to when any one 
of the beads reaches a switch position. A periodic solution for 
a system of A'^ particles will return to the initial configuration 
(both particle and trap positions) after 2N switches. 



where x{t) is the n component vector of the (local) posi- 
tion of the particles, s is the n-component vector of the 
(local) positions of the minima of potential proper for 
such time interval (i.e. s depends on the x{t) and on the 
previous history of the system), C„ is a real symmetric 
circulant matrix of order n. The first row of C„ (which 
defines the entire matrix) is: 



C„ = 



cos 27r/n -I- cos^ 7r/n cos 47r/n + cos^ 27r/n 



sin7r/n ' sin27r/ri 

cos 2(n — l)n/n + cos^(n — l)7r/n \ 
' sin(n — l)7r/n j 

III. THE DETERMINISTIC DYNAMICS 



(3) 



We now study the deterministic system described by 
Eq.([2]). The hydrodynamic coupling affects the over- 
damped oscillations, and determines the "phase differ- 
ences" between particles. Most remarkably, it leads to 
collective synchronized motions, which as we will now 
show are best understood by considering the eigenstates 
of the Oseen coupling matrix C„. 

Let us call Xj, j = 1, ■ ■ ■ ,n, the eigenvalues of the ma- 
trix C„ arranged in increasing values, and e^-'^ the corre- 
sponding normalized eigenvectors, C„ e*--'-' — Xj e^^\ Let 
us call hj(t) = [e^-^\x{t)) the projection of the position 
vector on the eigenvectors, that is the normal modes of 
the linear system. Projecting the system of Eq. [2] onto 
each of the eigenvectors gives: 



= 0, (4) 



and one easily obtains the uncoupled equations of motion 
for each normal mode: 



,(t)-(e<^)-F)+r,|/.,(t)=0 



with T., = 



3q \ 



(5) 



One sees that each normal mode hj{t) is related to 
a single exponential relaxation time Tj, and its linear 
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differential equation is homogeneous, if and only if the 
term (e^^^ • s) vanishes. In general, this is not the case. 

A solution of the dynamical system can be obtained 
with the following simple strategy (Fig[2|: Given initial 
positions x{ta) and initial configuration of potential s^°\ 
the evolution is trivially evaluated up to the time ti 
of the first hit. At this time one entry of the vector 
s changes sign, and the uncoupled system of Eq. [5] is 
again evaluated with the new constant vector s^^' up to 
the hit at time ^2- A basic role is then played by the 
continuity of each mode at the time of the hits: the fi- 
nal value before the hit becomes the new initial condition. 

The rotational discrete symmetry of the problem al- 
lows a detailed analytic study of the coupling matrix C„. 
The most relevant features are the following: 

(1) for every n being even integer, the lowest eigen- 
value Ai is a singlet, and its corresponding eigenvector 
ise<i) = ^(l,-l,l,... ,-1); 

(2) for every n being odd integer n > 5, the lowest eigen- 
value Ai is a doublet. The two eigenvectors may be writ- 
ten as 



^ f 1, cos(n — l)TT/n, cos 2{n — l)n/n, ■ ■ ■ 
cos(n — 1)^11 /nj , 

= —= f 0, sin(n — l)7r/n, sin 2(n — l)7r/n, • • • 
sin(n — Vf'TX jn ] ; 



(3) for every n, the matrix C„ has the eigenvector 
e — (1, 1, 1, • • • , 1). Only for n = 3 this constant 
eigenvector corresponds to the lowest eigenvalue Ai; 

(4) Since tr[C„] = 0, the lowest eigenvalue Ai < 0, and 
therefore the relaxation time ri > tq for every n. 



IV. THE ANALYTIC ASYMPTOTIC 
SOLUTIONS 

The dynamics of the system depends on three combi- 
nations of parameters: tq ~ 7o/k sets the unit of time 
for the evolution between hits; 3a/(8i?) measures the 
strength of the hydrodynamic interaction, and is limited 
by the bound d = 2_Rsin7r/n 3> 2a. Finally ^/A, with 
< ^/A < 1/2, fixes the amplitude of the oscillations. 
The periodic solutions are here obtained for arbitrary val- 
ues of the parameters (within the physical bounds given 
above). The plots for three and four particles, in Fig- 

m 



ures 
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^ ^ ^ have il\ = 1/4 and 3a/(8i?) = ^3/5. The 
plots for five particles, in Figures [7]and|8]have ,f/A = 1/4 
and 3a/(8i?) — 0.2. These values of 3a/(8i?) are much 
greater than sensible experimental values [3l |4j , but are 
chosen here to spread the relaxation times Tj, so empha- 
sizing the differences between periodic solutions. 



A. Case of n = 3 

For n — Z the lowest eigenvalue of the coupling matrix 
C3 is Ai = — a singlet. The next eigenvalue is a 

doublet A2 = A3 = The corresponding eigenvectors 
may be written: 




1 



\/2 \ _i 



e<3) 



1 

71 



(6) 



We find the following periodic solutions for the system, 
here listed in order of decreasing duration of the periods: 

(a) the three particles are synchronized in-phase; 

(b) a pair of equivalent solutions, where the trajectories 
are phase- locked (more precisely, time-shifted). 

Let us consider first case (a), with the configuration of 
particles in-phase^ i.e. such that xi{t) — X2(t) — 2:3 (i). 
The vector of the positions of the potential is 



s = ± 



The lowest normal mode is hi(t) = 
{xi{t) + X2{t) + X3{t)) ^/3xi{t) and it solves 

the equation: 

/ii(t)-(e<i)-s)+ri^/ii(i) = 0, 

rA , \/3a\ 




where ( e*-^ 



±y/S— , and ri = tq I 1 



8R 



Since (e^-'^ • s) = for j 7^ 1, the modes h2{t), h^{t) 
solve a homogeneous equation and the choice h2{t) = 
h^{t) = is consistent. As the three particles hit si- 
multaneously a geometric switch border, all of s — >■ —s. 
Therefore by continuity one obtains the opposite motion 
in the following time interval. The three particles per- 
form in-phase oscillations, satisfying the same equation 
as if they were un-coupled, as depicted on Fig. [IJb), but 
with the longer period 

— — log . 

2ti ^ ? 



Let us now obtain the phase-locked solutions, case (b) 
above. The symmetry of the problem, and results of nu- 
merical simulations, suggest to search for solutions where 
the time interval between hits is fixed, say A. Assum- 
ing an initial time to Sit one hit, the continuity of the 
trajectories requires: 



hj{to - 



rA) = 



)+ hj{to + ir-l)A)-{e 



gir)' 



with r = 1, 2, 3. 



(7) 
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0.2 0.3 0.4 



FIG. 4. (Color online) The period of the phase-locking solu- 
tion, dashed line, and the period of the solution with the 3 
oscillators in phase, in solid line, are plotted versus ^/A. The 
periods, in units of ro, were evaluated with 3a/{SR) = \/3/5, 
that is tq/ti — 0.8. 




The initial conditions of this periodic solution are then 



FIG. 3. (Color online) The fastest (shortest period) solution 
for n = 3, is phase-locked, (a) The trajectories of the particles 
1 (dash-red), 2 (solid-green), 3 (dash/dot-blue) are plotted 
versus time, for one period, (b) The normal modes hi{t) 
in red, /i2(i) in green, /i3(t) in blue are plotted versus time. 
Time is expressed in units of tq. This solution is not stable, 
as can be shown by numerical simulation, see section |IVE[ 
The other solution for n = 3 is the motion of all particles 
relaxing and switching in phase: it has longer period than 
the one illustrated here. 



By assuming that the system has a period T = 6 A, one 
obtains: 

We choose to when a;i(io) = —A/2 -|- ^, the particle 2 
is also moving toward A/2 and particle 3 is moving in 
the opposite direction. That is we consider the sequence 
{s(j)}, with j = I, • • • , 6 , s<^'+3) = -s<J) , .s<i) = s<^): 





Finally by inserting the initial condition 
xi{ta) = ^hiito) - ^hsito) = -f + we find 
the equation that fixes A in terms of the parameters 
of the problem. This is described in Appendix. Fig[3^ 
shows the trajectories of the particles 1, 2 and 3 as 
a function of time, for one period T = 6A. Time 
is expressed in units of tq. Figj3jD plots the normal 
modes hi{t), h2{t) and h^{t)^ again for one period. The 
oscillations of the pair of normal modes corresponding 
to the degenerate eigenvalue have an amplitude much 
larger than the oscillations of the normal mode hi{t). 
We verified that there exists a completely analogous 
phase locked synchronized solution where the role of 
particles 2 and 3 are exchanged. 

In Fig|4| we compare the period of the phase-locking 
solution, with the period of the solution with the 3 oscil- 
lators in phase, as a function of ^/A. One sees that the 
solution where the three particles are synchronized in- 
phase always has a longer period than the period of the 
pair of phase-locked solutions, independently of the geo- 
metric conditions. This is also the solution which appears 
in the experiment and in the numerical simulations (4j. 
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B. 



Case of n = 4 



For n = 4 the lowest eigenvalue of the matrix C4 is 
a singlet Ai = —1 — \/2, the next is also a singlet A2 = 
— 1 + v^, and the next two are a doublet A3 A4 = 1. 
The eigenvectors may be written: 




(10) 



The periodic solutions of the system in order of 
decreasing periods are : 

(a) adjacent particles in anti-phase configuration; 

(b) all particles synchronized in-phase; 

(c) the pair (1,2) are anti-phase with the pair (3,4). 
Also an equivalent solution where the pair (1,4) are 
anti-phase with the pair (2,3); 

(d) a pair of phase-locked solutions, with the same period. 

In the anti-phase configuration (a) , the vector s of the 
positions of the potentials is proportional to the eigenvec- 
tor e^^) . Then, all the normal modes hj (t) with j = 2,3,4 
obey homogeneous differential equations and may consis- 
tently vanish at all times. The normal mode hi(t) gen- 
erates the following solution: 



(a) 

0.2 

x(t) 




/ 1 




(b) 

0.2 

x(t) 

0.1 



0.1 
0.2 




1.5 t 




FIG. 5. (Color online) (a) The "slowest" solution (longest 
period) for the case n — 4 has nearest neighbors in antiphase 
with each other (hence next-nearest neighbors are in-phase). 
The trajectories of the pair of particles 1 and 3 in red-dash, 
and the other pair in green-solid are plotted versus time, in 
units To. (b) There are two equally "fastest" solutions for 
n — 4. One has nearest-neighbor pairs in phase, and the other 
pairs in antiphase with each other. This would look like panel 
(a), but with particles 1,2 in red, and 3,4 in blue. The other 
solution is shown in (b), and has trajectories of the particles 1 
to 4, respectively in red-dash, green-solid, blue-dash/dot and 
cyan-thick are plotted versus time, in units tq. 



s = ± 



A -1 



±A, 



-1 



/ii(i)= (e<i).f(i)) -2a;i(<), 

Xi{t) = -X2it) = X3{t) ^ -X4{t), 

/ii(i)- +Ti|/ii(t)=0, 
(l + \/2)3aV^ 



with Ti = To 1 



8R 



Here the four particles perform anti-phase oscillations, 
satisfying the same equation as un-coupled oscillators, 
but with Ti replacing tq. The period of this anti-phase 
solution is 



Ti = 2ti log 



The trajectories for this case are shown in Fig. [5^. 

The next solution, case (b) above, is the in-phase syn- 
chronization and may still be written in terms of just one 



normal mode, h2{t). Its relaxation time and period are 



T2 = To 



8R 



T2=2T2log^^<ri. 

In the next solution, case (c), the couple xi{t) — X2{t) 
are in anti-phase with the couple x^it) — 2:4 (t). The 
vector s of the position of the potential is proportional 
to the sum e^'^' -|- e^*^^ . Then we may choose hi (t) = , 
h2{t) = 0, and: 



s = ± 




V2' 



Xi{t) = X2{t) = ~X3{t) = ~X4{t) = --j=hi{t), 

halt) = hi{t), 



T3 = To 1 



3a 
8^ 



Ts = 2t3 log ^^r^ <T2. 
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The trajectories perform oscillations as depicted in eigenvectors may be written 
Fig. [5^, but the relaxation time is replaced by T3 and the 
pair of particles 1 and 2 in red, the other pair in blue. A 
completely equivalent solution has the pair of particles 
xi{t) = X4{t) in anti-phase with the pair X2{t) — 0:3 (i). 



Finally there are phase-locked solutions, case (d) 
above, also with short periods T, where X4{t) = xi{t — 
3r/4), X3{t) = xi{t - r/2), X2it) = xi{t - T/4). A 
completely equivalent periodic solution has the reverse 
order Xi{t) ^ xi{t + 3T/4), xsit) = xiit + T/2, X2{t) = 
xi(t + T/4). 

Since the sequence of the positions of the potentials is 
s = I (1, -1, -1, 1) -s -)• s -)■ -s -)■..., we find 
hx[t) = , /i2(t) = 0, and 



with r3 = 



A (1 



-A/T3' 



V2 l + e-2A/r3 
71 1 + e-2A/r3 ' 



1 



3a 
8i? 



(11) 



The period is: 




T = 4A = 2t3 log^-^ =r3. 



where ^ = A 



1 + e-2A/r3 ' 



(12) 



same period as case (c), and the motion is given by: 



xi{t) 



'V2 



hi{t), X2{t) 



\/2 



The trajectories for this are shown in Fig. [5]3, for the 
particles 1 to 4. 



The periodic solutions of the system in order of 
decreasing periods are: 

(a) a pair of phase-locked solutions, one with 
xi{t) = X3{t + 2T/10) = X5{t + 4T/10) = X2{t + 6T/10) = 
X4{t + 8T/10) and the analogous solution with opposite 
time-shifts; 

(b) a pair of phase-locked solutions, one with 
xi{t) X2{t + 2T/10) = .T3(i + 4T/10) = a;4(< + 6r/10) = 
x^it -\- 8T/10) and the analogous solution with opposite 
time-shifts; 

(c) a synchronized in-phase solution Xj{t) ~ x{t) , 
J = l,2,--- ,5. 

The periods of these solutions are shown in Fig. [6] and 
evaluated in the Appendix. 



C. Case of n = 5 



The odd-numbered systems, except rt = 3 discussed 
above, have general properties. For n = 5 the lowest 
eigenvalue of the matrix C5 is a doublet Ai — A2, the 
next is a singlet A3, followed by a doublet A4 = A5. The 



The phase-locked solution (a) is quite relevant since 
for every odd number of particles greater than 3, the 
system synchronizes into a periodic solution of this type. 
At any time, three particles are moving in one direction 
and 2 are moving in the opposite direction. We call A 
the uniform time interval between any two consecutive 
hits. The period T — 10 A. Take, for example, that 
at time to the particles 1,3,5 are moving towards the 
positive direction, with Xiito) < X3{to) < x^itQ). Then 
the sequence of the vectors {s^-'-'}, with j = 1, • • • , 10 , 
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T 

3.5 - 



3a 
8~R 



FIG. 6. (Color online) The periods for the phase-shifted solu- 
tion (red-dash) with adjacent trajectories {1, 3, 5, 2, 4}, for the 
phase-shifted solution (green-solid) with adjacent trajectories 
{1, 2, 3, 4, 5}, and for the in-phase solution (blue-dash/dot), in 
units of tq, are plotted versus the hydrodynamic coupling 
3a/(8-R). These periods are evaluated in the Appendix. 



Sin 



_s(j+5) 



-1 
1 

-1 
V 1 J 



1 

-1 



( w 



-1 
1 
1 

V-i/ 



(\\ (-}\ 



\-iJ 



A 



V-i/ 



(13) 

No eigenvector e^^ is orthogonal to all of the eigenvec- 
tors {s^-'-*} of this sequence, therefore all the five normal 
modes hj{t) contribute to the solution. 

The continuity and periodicity conditions fix the initial 
conditions: 



,(^o)( 



1 



-5A/T, 



{5-r)A/rj ^ 



= -(l-e^/-^)^(e<^).s<'-)j e 

r=l 

I 1 - e-5A/r, 

-l + 2e-2A/r, „g-5A/r, 

1 - 2e~^'^/'^J -I- e^^^l'^i 

and as usual r, = 
One last equation determines the time interval A. If 

A 



(14) 



we use x\(ti^ = — 2 + it is 



(/ii(to) + /i4(io)) 



/i3(to). (15) 



The trajectories are shown in Fig. [7^, for the particles 
1 to 5 (the coupling is chosen at 3a/(8i?) = 0.2, and 




FIG. 7. (Color online) This phase locked solution for n = 5, 
with the nearest neighbors almost in anti-phase, is the so- 
lution with the longest period, and is the stable state with 
attractive potentials, (a) The trajectories of the particles 1 
to 5, respectively in red, green, blue, red-dash and green- 
dash are plotted versus time, in units ro. To draw the figure, 
the coupling was chosen 3a/(8i?) = 0.2, then the period of 
the solution is T ~ 3.8547ro. (b) The normal modes /ij (i) , 
j = 1,2, •••5, respectively in red, green, blue, red-dash and 
green-dash are plotted versus time, in units ro. 



then the period of the solution is T ~ 3.8547ro). 



Fig. [Tja shows the complex behaviour of normal modes 
h.j{t) , j = 1,2, •••5. The doublet of normal modes 1 
and 2 present oscillations with an amplitude much larger 
than the oscillations of the doublet of normal modes 4 
and 5. The singlet normal mode 3, related to the center 
of mass X] j=i even smaller oscillations. 



Another pair of phase-locked solutions, case (b), 
with shorter period has the trajectories of particles 
{1,2,3,4,5} and the equivalent solution with the ar- 
rangement of trajectories {1, 5, 4, 3, 2}. Here, like in case 
(a), at any time there are three particles moving in one 
direction and 2 moving in the opposite direction. We call 
A the uniform time interval between two hits. The pe- 
riod r = 10 A. Taking for example that at time io the 
particles 1,4, 5 are moving toward the positive direction, 
with xi(io) < 2:5 (to) < 2:4(^0), then the sequence of the 
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vectors {s^-?'^}, with j = 1, • • • , 10 , s^^'^ = ~s^i+^) is 



2^1 
A 



/ \ \ 



V 1 / 



V 1 / 



V 1 / 



-1 
-1 



1 
1 

-1 

V-i/ 



1 
1 

-1 

\-lJ 



(16) 



As before, the continuity and periodicity conditions fix 
the initial conditions: 



and as usual r,, = 



/ 1 ~ 

1 - 2e"^'^/^J' + e-^'^/^^ 
1 - 2e-^/^^ + e-5A/rj 
-1 + 2e"''^/^J - e^^'^/^J' 

To 



1 



3a \ 



(17) 



By using the same initial condition, a;i(io) = — 2 + 
the equation that determines A still has the form of 
Eq. (15 1. The trajectories are shown in Fig. [sk, for 



the particles 1 to 5. In the figure, the coupling is 
3a/(8i?) = 0.2, and then the period of this solution is 
T - 1.7407ro. 

Fig. [sJd shows the normal modes hj{t) ,j — 1, 2, • • • 5. 
The doublet of normal modes 4 and 5, which have the 
shortest relaxation time, oscillate with much larger 
amplitudes than either the doublet of normal modes 1 
and 2 or the singlet normal mode 3. 

The synchronized in-phase solution, case (c), for the 
five particles is very simple, since only the mode h^lt) 
contributes, whilst all the hj{t) = for j 7^ 3. The 
relaxation time and period are 



T3 = To 1 



3a 
8^' 



A. 



= 2t3 log 



(18) 



This is the "fastest" periodic solution. 



D. Case of higher n 

As the number of particles of the system increases, 
the analytic periodic solutions of the system become 
more and more complex. Yet some few simple features 




FIG. 8. (Color online) This phase-locked solution has sorter 
period than the one depicted in Fig[7| (a) The trajectories of 
the particles 1 to 5, respectively in red, green, blue, red-dash 
and green-dash are plotted versus time, in units tq. To draw 
the figure, the coupling 3a/ (87?) = 0.2, then the period of the 
solution is T ~ 1.7407to. (b) The the normal modes hj{t) 
,j = 1, 2, ■ • ■ 5 , respectively in red, green, blue, red-dash and 
green-dash are plotted versus time, in units tq. To draw the 
figure, the coupling 3a/{8R) — 0.2. 



are generic. 

We already mentioned that for every particle number 
n being an even integer, the lowest eigenvalue Ai is a 



singlet, its eigenvector is e^^^ = (1, — 1, 1, • • • ,— 1). 
This allows a periodic solution where adjacent particles 
are in anti-phase, just as the first solution (a) described 
for the system of four particles. Since hi(t) is the unique 
normal mode non- vanishing, the period of the anti-phase 
solution is 



Ti =2ri log^-^, Ti =To ( 1 



n/2-1 

Ai = (-l)"/2+i_2 (-1) 



3a 
8^ 



Ai 



where 



.cos^-Hcos^i^ 



sm 



JIT 



that is, 



in particular: 



Ai 2.414 (n = 4); 



Ai - -4.577 (n = 6); 



Ai ~ -6.528 (n = 8). 
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This anti-phase solution is the periodic solution with the 
longest period, for every even n. 

We also mentioned that for every n being odd integer 
n > 5, the lowest eigenvalue Ai is a doublet. The 
periodic solution with the longest period is analogous to 
the solution (a) for the system of five particles. That 
is a pair of equivalent solutions where a sequence of 
trajectories are time shifted by a uniform time interval. 
For n = 7 the pair of sequences is {1,3,5,7,2,4,6} 
and {1,6,4,2,7,5,3}. One may notice that if n » 1, 
adjacent particles are time shifted by almost half period, 
that is are approximately in anti-phase. Thus the 
distinction between even and odd number n decreases, 
as n increases. 



x{t) 




FIG. 9. (Color online) Trajectory of a driven oscillator with 
repulstve driving force, uncoupled to other oscillators. Note 
the increase in velocity between switch positions, characteris- 
tic of the motion with potentials of negative stiffness. 



E. Simulations 

To verify which of the analytical solutions is the dy- 
namical steady state, we resorted to numerical integra- 
tion. Two different codes have been developed. 

Firstly we solved the coupled system of Eq[2] which 
corresponds to the tangential motion, using the solver 
odell3 of matlab, which implements the Adams- 
Bashforth-Moulton multi-step (i.e. timestep optimizing) 
method. We used the event handling facility integrated 
into all matlab 's ode solvers to intercept the instant of 
time at which any one of the particles hits the switch 
position. As a stabilizing measure, we imposed that at 
a "hit" all particles which are in the range of Ss from 
their switches are considered to have hit their target and 
consequently their centers of force switch their positions; 
several values in the interval 10~^ X < Ss < 10~^ A have 
been used with no substantial difference on the results. 

Secondly, a different code described in [3l HI [14] with 
fixed timestep, is also used, enabling testing in the ab- 
sence of noise but also the robustness at finite tempera- 
ture. As well as exploring noise, this Brownian Dynamics 
(BD) code solves a model which is very close to the ex- 
perimental condition [3] : the colloidal particles are free to 
move in a two dimensional environment, and the tangen- 
tial and normal forces can be set independently at differ- 
ent stiffnesses. Finally, experiments typically work at a 
finite sampling rate, and this can be mimicked closely in 
the fixed-timestep BD code as was discussed in previous 
study j3]. 

The numerical solution (both the simulation ap- 
proaches described above) show the remarkable result 
that for every particle number n and every initial condi- 
tion, the system evolves to a synchronized configuration 
corresponding to the periodic solution with the longest 
period (for the case of positive trap stiffness k > dis- 
cussed up to here). If such a configuration is a pair of 
equivalent solutions, as it happens for odd n > 5, the 
synchronization occurs on either one of the pair. These 
findings also agree with the experiment in [1]. 



Another result of the extensive numerical solutions is 
that we never observe different solutions from the cases 
that were investigated analytically (in principle, we could 
not have ruled out the presence of much more complex 
and less symmetric solutions). 



V. REPULSIVE POTENTIAL 

It is interesting to consider a harmonic repulsive po- 
tential, K < 0, acting on the colloidal spheres. This may 
be realized experimentally by tailoring an appropriate 
potential landscape with time-sharing or holographical 
optical traps. Then the synchronization of the system 
of oscillating particles due to the hydrodynamical 
interaction still occurs, but on a periodic configuration 
different from the one observed in the attractive case. 
With minimal changes from the previous analysis we 
can describe the periodic solutions occurring in the case 
of repulsive harmonic potentials. Again the numerical 
simulations of this case confirm that the systems con- 
verge onto one of the analytical solutions, and in this 
case it is one with short period. 

With negative stiffness, if we ignore the hydrody- 
namic coupling, a particle moves toward the pos- 
itive X axis according to the differential equation 
7o i{t) + K (x{t) + f) — f{x) = 0, where now k < is the 
stiffness of the repulsive harmonic force, tq = 7o/|k|, / 
is a stochastic force due to thermal agitation of the fiuid, 
which will be ignored in the analysis of the deterministic 
system. In this case, the motion is exponentially accel- 
erated until the particle reaches the position ^ — At 
that time the center of the repulsive harmonic potential 
moves from the position — to ^. The particle inverts 
its motion until it reaches the position — f + C when 
the potential jumps again. The resulting oscillatory 
motion for a single particle (or one in a non-coupled sys- 
tem), neglecting the stochastic force, is depicted in Fig.|9] 

The period Tq of these uncoupled oscillations is 
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To = 2to log . The oscillatory motion is similar to 
the attractive potential case, shown in Fig. [T] except for 
the concavity. 

One may repeat the analysis of the attractive case, and 
obtain the linear deterministic system, valid in a time 
interval between two hits: 

/ + ||c„) m Topit) = 0, TO = 

(19) 

where x{t) is the n component vector of the (local) posi- 
tion of the particles, s is the n component vector of the 
(local) position of the minima of the repulsive potential 
proper for such time interval (opposite of the previous 
case), C'n is the same real symmetric circulant matrix 
of order n, given in eq.dsl). We again define the normal 
modes hj{t) = (e^-'^ x(i)j and their decoupled equations: 



0, 



•3 



1 



3a \ 

d>R^3 



(20) 

They differ from the previous equations ^ only for 
the sign of the time evolution. It follows that the 
set of periodic solutions, described in the previous 
analysis for the attractive traps, still exist in the 
present case of repulsive traps, with the same periods. 
Only the shape of the trajectories is affected by re- 
placing diverging exponential in place of converging ones. 



Fig. [TOfi and Fig. [TOp are examples of this point, 
when compared with the Fig. [3^ and Fig. [Sja. We 
plotted the phase-locked solution for three particles in 
the repulsive case. Fig. [TO^ shows the trajectories of 
the three particles as function of time, for one period 
T = 6A. Fig. 10 3 plots the normal modes hi{t), h2{t) 
and h^{t)^ again for one period. The oscillations of the 
pair of normal modes corresponding to the degenerate 
eigenvalue have an amplitude much bigger then the 
oscillations of the normal mode hi{t). Remarkably, the 
simulations show that in the case of repulsive traps it is 
this periodic solution, the one with the shortest period, 
which is the stable asymptotic synchronization for the 
system. 

With systems of higher n, we have observed in the 
numerical simulations that again on switching from at- 
tractive to repulsive traps the stable solution ceases to be 
that with longest period. Instead, the system converges 
onto one of the solutions with shortest period. However, 
as is clear by observing the sequence of solution peri- 
ods T calculated for the cases (n = 4 or n = 5), while 
there is a clearly separate largest T, there are various 
solutions with small T. We have tested numerical solu- 
tions in the physically realistic parameter space, using 
the two numerical approaches outlined in section |IVE| 
the ODE solver with the "strong" tangential motion con- 
straint converges to the shortest period solution, whereas 




FIG. 10. (Color online) (a) The trajectories for the the phase- 
locked solution for three particles in the repulsive case are 
plotted versus time, in units tq. (b) The normal modes hi(t) 
in red-dash, h2(t) in green-solid, h^lt) in blue dash/dot are 
plotted versus time, in units of tq. 



the Brownian Dynamics simulation with the weaker "har- 
monic" constraint converges to one of the shorter period 
solutions, but not always the shortest one. This will be 
explored further in future work; our current understand- 
ing is that the system with repulsive potentials converges 
to one of the solutions with shorter periods, and that 
when there are various periods grouped close together 
the details of the system play a strong role in selecting 
between these. 



VI. CONCLUSION 

The model studied in this paper has several unusual 
features. Most studies on synchronization show that it is 
a threshold process: it occurs when the coupling among 
the oscillators is strong enough or when increasing the 
number of oscillators, a large fraction synchronize to 
a common frequency |15j . In contrast, in the present 
model, where we studied the deterministic system 
without stochastic forces, we find no threshold. Synchro- 
nization occurs for any strength of the hydrodynamic 
coupling, that is any value of a/i?, and any number of 
oscillators. 

Most models on synchronization are described by 
non-linear differential equations where a very limited 
progress is achieved analytically, whereas the present 
model is piece-wise linear and all periodic solutions may 
be derived. 

Numerical solutions in the case of attractive harmonic 
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potentials clearly show the synchronization to the peri- 
odic solution with the longest period, for any number of 
particles. In the case of repulsive harmonic potentials, 
numerical simulations show the synchronization to the 
the pair of phase locked solutions in the case of three 
particles. We have not performed a stability analysis 
about the periodic solutions, but there is no reason 
to expect disagreement with the results of numerical 
simulations. This change of the asymptotic solution 
with the change of curvature of the potential had 
already been shown numerically for a linear chain in |12| . 
The analytical framework developed here adequately 
describes experiments on systems of colloidal particles 
with moving optical traps [4]. 



VII. APPENDIX 

In the case of three particles, with attractive harmonic 
traps, the initial conditions for the normal modes, for 
periodic phase-locked solutions are given in the cqs.([9]). 
We supplement them with an initial condition: 

12 A 
and obtain an equation that determines A: 

+ ^ (1 - e-'^/^^) (e-2A/r2 + 2e-^/^^ + l) + e^'^^'^) 
= 0-e)(l + e-3A/n)(i + ,-3A/..y (21) 



Furthermore, tq/ti +2to/t2 = 3, since Ai -I-2A2 = 0. Let 
us define: 



Pi 



To 



< 1, P2 



To 
T2 



>1, Pi+ 2p2 = 3, 



then Eq.(|21|) may be rewritten in the simpler form: 

„2pi _ „3p 



A r 
3 



1 + xP^ (1 -I- xP^ 



2x^P' 



{-2-xP' +x^P'-3x^P'] 



„3pi 



1-X 



P2 



= 0. 



(22) 



If we choose the coupling strength 3a/(8i?) = a/3/5, we 
have pi — 0.8 and p2 — 1.1. The further behaviour of A 
as function of A/^ is depicted in Fig. [4] 
If we further choose ^ = A/4, then eq.pTj) implies x ~ 
0.70701565, that is: 



A 

To 



log x - 0.346702, T = 6A = ~6tq log x . 



For the case of five particles with attractive traps. 
Fig. [6] shows the period T in units of tq for the phase- 
shifted solution with adjacent trajectories {1,3,5,2,4}, 
the period of the phase-shifted solution with adjacent 
trajectories {1,2,3,4,5}, and the in-phase solution, ver- 
sus the hydrodynamic coupling 3a/ (Si?), while we fixed 
A = 1, ^ = 0.25. In the limit of vanishing coupling, the 
three curves converge to T = 2tq log(A— ^)/^ ~ 2.1972 tq. 
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